Real space finite difference method for conductance calculations 
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We present a general method for calculating coherent electronic transport in quantum wires and 
tunnel junctions. It is based upon a real space high order finite difference representation of the single 
particle Hamiltonian and wave functions. Landauer's formula is used to express the conductance 
as a scattering problem. Dividing space into a scattering region and left and right ideal electrode 
regions, this problem is solved by wave function matching (WFM) in the boundary zones connecting 
these regions. The method is tested on a model tunnel junction and applied to sodium atomic wires. 
In particular, we show that using a high order finite difference approximation of the kinetic energy 
operator leads to a high accuracy at moderate computational costs. 

PACS numbers: 73.63.-b, 73.40.-c, 71.15.-m, 85.35.-p 
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I. INTRODUCTION 

The progress in experimental control on the nanometer 
scale has enabled studies of electronic transport in quan- 
tum wires of atomic dimensions^ The transport proper- 
ties of such systems have to be understood on the basis 
of their atomic structure. This notion has generated a 
large effort in recent years to calculate the conductance 
of quantum wires from first principles. Several different 
approaches have been formulated, which have a common 
basis in the Landauer-Buttiker approach to express the 
conductance of a coherent system in terms of a quan- 
tum mechanical scattering problem^ In such calculations 
the quantum wire consists of a scattering region of fi- 
nite size, sandwiched between two semi-infinite leads that 
arc considered to be ideal ballistic wires. Semi-empirical 
tight-binding models have been exploited to solve this 
problemA^i& Aiming at a better description of the elec- 
tronic structure, several current approaches rely upon 
density functional theory (DFT). 

The main differences between these approaches lie in 
the approximations that are used to describe the atomic 
structure of the leads and in the techniques that are used 
to solve the scattering problem. In pioneering work, jel- 
lium (i.e. free electron) electrodes have been used to de- 
scribe the leads and the scattering wave functions have 
been obtained by a transfer matrix method 7 - or by solv- 
ing the Lippman-Schwinger equation*^ A transfer ma- 
trix method has also been used taking into account the 
full atomic structure of the leads at the DFT level>i2iii 
Alternatively, the conductance can be calculated using 
a Green function approach without calculating the scat- 
tering wave functions explicitly^ Several implementa- 
tions of this approach have been formulated that use a 
localized basis set to form a representation of the scat- 
tering problem. These implementations mainly differ in 
the kind of basis set used, e.g. Gaussian or numerical 
atomic orbitals, or wavelets ji&iii&ii&iLiLiili An embed- 
ded Green function approach has been applied using a 
delocalized basis set of augmented plane waves' 

In this paper we present a technique for solving the 



scattering problem of a quantum wire without the use 
of a basis set. Instead, potentials and wave functions 
are represented on a uniform real space grid and differ- 
ential operators are approximated by a finite difference 
approximation (FDA). Previous implementations of this 
idea have used a simple first order FDAi 21 ' 22 : 2 ?' 24 In that 
case the grid has to be relatively fine in order to ob- 
tain sufficiently converged results. This hinders the ap- 
plication to large systems because of the computational 
costs involved in using fine grids. However, in ground 
state (DFT) electronic structure calculations high order 
FDA's have been shown to markedly increase the effi- 
ciency of real space grid techniques by enabling the use 
of coarse grids. 25 ' 26 27 In this paper we demonstrate that 
high order FDA's make it possible to solve the scattering 
problem much more efficiently. 

The method we propose for calculating the conduc- 
tance of a quantum wire is based upon wave function 
matching (WFM) in the boundary zones connecting the 
leads and the scattering region. 28 Unlike transfer matrix 
methods, however, it does not require the explicit calcu- 
lation of wave functions in the scattering regioniifiiii 
It does not require the explicit calculation of Green 
functions either) 12 ! 1 ?! 14 ! 15 : 16 ! 1 !: 18 ! 1 ?] 2 ? which enables us to 
solve the scattering problem at real, instead of complex, 
energies^ Our method can be classified as an O(N) tech- 
nique, since the computing costs are determined by the 
size of the scattering region with which they scale lin- 
early. A related technique that uses a linearized muffin 
tin orbital basis set, has been applied to calculate the 
electronic transport in layered magnetic materials 
Although the formalism presented here can be extended 
to the non-equilibrium situation, we consider in this pa- 
per the linear response regime only. 

This paper is organized as follows. In Sec. [H] the 
main ingredients of our computational method are ex- 
plained, where the computational details can be found in 
appendix ^ The accuracy and convergence properties 
of the method are verified on model tunnel junctions in 
Sec. IIII Al The application to a more complex system, 
which consists of a sodium atomic wire, is discussed in 
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Sec. IIII Bl A summary is given in Sec. 



II. COMPUTATIONAL METHOD 

Within the Landauer-Biittiker approach the conduc- 
tance G of a quantum wire is expressed in terms of the 
total transmission T(E) 



(1) 



assuming spin degeneracy^ T(E) can be obtained by 
solving the quantum mechanical scattering problem at 
the fixed energy E. Eq. is valid in the linear re- 
sponse regime, where T{E) needs to be evaluated at the 
Fermi energy E = Ef- Our quantum wire is defined as a 
system consisting of a finite scattering region that is con- 
nected left and right to semi-infinite leads. The latter 
are supposed to be 'ideal' wires, which can be described 
by a periodic potential along the wire direction. In the 
scattering region the potential can have any shape. We 
consider two cases that can be treated by essentially the 
same technique. In the first case the system has a finite 
cross-section perpendicular to the wire direction, whereas 
in the second case the system is periodic perpendicular 
to the wire. The latter case also covers planar interfaces 
and tunnel junctions. 

In order to solve the scattering problem we generalize a 
method formulated by Ando.— Here one basically solves 
a single particle Schrodinger equation directly at a fixed 
energy E in two steps. In the first step one obtains the 
modes of the ideal leads. Subsequently the wave func- 
tions for the scattering region are constructed such, that 
they are properly matched to the solutions in the leads. 
We use a real space finite difference method to represent 
the Schrodinger equation. In the following three subsec- 
tions we will introduce this representation and discuss 
the steps required to solve the scattering problem. 



A. Finite difference approximation 

We start from a single particle equation of the general 
form 



E - V(r) 



2m 



*(r) = 0, 



(2) 



which represents the Schrodinger equation of a single par- 
ticle in a potential V. Alternatively, within the DFT 
scheme it represents the Kohn-Sham equation with V 
the total effective potential. We put the wave function 
^ and the potential V on a equidistant grid in real space 
r = (xj,y k ,zi), where Xj = x +jh x ,y k = ya + kh v , z t = 
zq + lh z and h Xl h y , h z are the grid spac ings in x, y and z 
directions, respectively. Following Refs. I25l26l we replace 
the kinetic energy operator in Eq. (J2J by a high order 
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FIG. 1: (Color online) (a) The system is divided into cells 
indicated by an index i. The cells have L ■ W y ■ W z grid 
points in the x, y and z directions, respectively, ^ffi is the 
supervector that contains the wave function values on all grid 
points in cell i. (b) Hi is the Hamilton matrix connecting 
grid points within cell i; the B-matrix connects grid points 
between neighboring cells and is independent of i. 



FDA. For the x part this gives 
1 N 

T2 c n^{x ]+ „,y k , 



d 2 ^(x J1 y k ,zi) 



dx 2 



z 



(3) 



n=-N 



with similar expressions for the y and z parts. Expres- 
sions for the coefficients c n for various values of N are 
tabulated in Ref. l2tj The simplest approximation in 
Eq. © (N = 1, where C\ = c_i = 1 and c = —2) reduces 
Eq. (J2J to the well-known simple finite difference repre- 
sentation of the Schrodinger equationii2i2ii22i2&2£ How- 
ever, we will demonstrate that the scattering problem 
can be solved much more efficiently using higher order 
FDA's with N = 4-6. 

In a FDA the Schrodinger equation of Eq. (J2J becomes 

JV 

- 0, (4) 



n=-N 



where V, &j,k,l is a short-hand notation for V, &(xj, y k , Zi) 
and tn V ' z = h 2 /2mh 2 , y z x c n . In order to make a connec- 
tion to Ando's formalism^ we divide the wire into cells 
of dimension a x x a y x a z . The direction of the wire is 
given by the x-axis. The number of grid points in a cell 
is L = a x /h x , W v — a y /h y , W z = a z /h z for the x, y and 
z directions respectively. We wish to distinguish between 
two different cases. In the first case the wire has a finite 
cross-section in the yz plane. In the second case the wire 
has an infinite cross-section, but it has a periodic poten- 
tial in the yz plane, i.e. V Jtk+Wy ,i = V jt k,i+w z = Vj,k,l- 
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In both cases the (unit) cell in the yz plane is described 
by W y x W z grid points. 

The values *&j,k,l where the indices j, k, I correspond to 
a single cell i are grouped into a supervector \Jv The idea 
is shown in Fig. ^ This supervector has the dimension 
N IS = L ■ W y ■ W z , which is the total number of real space 
grid points in a cell. If we let i denote the position of the 
cell along the wire then Eq. (0} can then be rewritten as 

(M-H i )* i + B* i _i+Bt* i+1 =0, (5) 

for i = — oo, . . . , oo. Here I is the -/V rs x _/V rs identity 
matrix. The matrix elements of the N IS x N IS matrices 
H, and B can be derived straightforwardly from Eq. 10}. 
The expressions are given in appendix IA II both for a 
wire that is finite and for a wire that is periodic in the 
yz plane. For the latter Hi = Hj(k|i), where k|| is a 
wave vector in the two-dimensional Brillouin zone. In 
the following this notation is suppressed. 

Eq. JSJ has the form of a nearest neighbor tight-binding 
equation, expressed in terms of vectors/matrices of di- 
mension iV rs . This form enables us to use Ando's tech- 
nique to solve the scattering problem^ Note however 
that the matrices B, B^ in Eq. Q are singular, see 
Eq. (|A3jl . which requires a generalization of this tech- 
nique. 

B. Ideal wire 

An ideal wire is defined by a potential that is periodic 
in the direction of the wire, i.e. V{xj + a x ,yk,zi) = 
V(xj,yk,zi) or V 3+L ,k,i = V 3 . k ,i- Since the potential is 
the same in each cell, the matrix H; = H in Eq. (J5J is 
independent of the cell position i. In a periodic system 
the vectors in subsequent cells are related by the Bloch 
condition 

A* 4 = * t+ i, (6) 

where A = e lkmax with k x real for propagating waves 
and complex for evanescent (growing or decaying) waves. 
Combining Eqs. (J5J and (JSJ one then obtains the follow- 
ing generalized eigenvalue problem 

[(V ?)-*(-?'!)](£ ) = »■ m 

Formally, the dimension of this problem is 2N rs . There 
are a number of trivial solutions, however, since B, B^ are 
singular matrices. In appcndix IA 2l it is shown how reduce 
the problem to its 2N ■ W y ■ W z non-trivial solutions. 

The non-trivial solutions of Eq. (JJJ can be divided into 
two classes. The first class comprises Bloch waves prop- 
agating to the right and evanescent waves decaying to 
the right; the corresponding eigenvalues are denoted by 
A(+). The second class comprises Bloch waves propagat- 
ing to the left or evanescent waves decaying to the left; 
the eigenvalues are denoted by A(— ). The eigenvalues of 



the propagating waves have |A(±)| = 1 and for evanes- 
cent waves |A(±)| ^ 1. The evanescent states come in 
pairs, since it is easy to show that for every solution A(+) 
there is a corresponding solution A(— ) = 1/A*(+). It can 
be shown that the propagating states also come in pairs, 
i.e. for every right propagating wave A(+) there is a left 
propagating wave A(— )Ji 

It makes sense to keep only those evanescent waves for 
which 1/5 < |A| < S, where 5 is a sufficiently large num- 
ber. States with |A| outside this interval are extremely 
fast decaying or growing. Such states are not impor- 
tant in matching an ideal wire to a scattering region. 
Typical of finite difference schemes there are also non- 
physical solutions to Eq. Q, which are related to so- 
called parasitic modesif&M These are easily recognized 
and discarded since their |A|'s are either extremely small 
or large and thus fall outside the selected interval. More- 
over, these |A|'s are very sensitive to the grid spacing and 
rapidly go to or oo if the grid spacing is decreased. 

After filtering out the physical and useful solutions 
Eq. Q) we end up with M pairs of solutions A m (±); m = 
1,...,M, where usually M <C N ts . We construct the 
normalized vectors u m (±) from the first N IS elements of 
the eigenvectors of Eq. (Q) and form the N IS x M matrices 

U(±) = ( Ui (±)---um(±)). (8) 

Choosing the cell i = as the origin, one then writes the 
general solution \I/o in this cell as a linear combination 
of these right- and left going modes 

*o = *o(+)+*oR, (9) 

where 

M 

* (±) = U(±)a(±) - ]T u m (±)a m (±), (10) 

with a(±) vectors of arbitrary coefficients of dimension 
M. 

Defining the M x M diagonal eigenvalue matrices by 

(A(±)) nm = S nm \ m (±), (11) 

and using the Bloch condition of Eq. the solution in 
the other unit cells then can be expressed in a compact 
form 

= U(+)A*(+)a(+) + U(-)A*(-)a(-). (12) 

In order to apply Ando's formalism^ it is advanta- 
geous to slightly rewrite this. We define the N IS x N IS 
matrices F(±) and F(±) by 

F(±)U(±) = U(±)A(±), (13) 
F(±)U(±) = U(±)A- 1 (±). (14) 

Note that F(±) ^ F _1 (±) since the N rs x M matrices 
U(±) are not square (typically M iV rs ). This presents 
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FIG. 2: (Color online) (a) Schematic representation of a quan- 
tum wire. The left (L) and right (R) leads are ideal wires that 
span the cells i = — oo, . . . , and i — S + ... , oo, respec- 
tively. The scattering region spans cells i = 1, . . . , S. (b) The 
reduced problem spans the cells i = 0, . . . , see Eq. (1 12 ( i I . 



no problem, however, and explicit expressions for the ma- 
trices F, F are given in appendix IA 31 They allow Eq. i|12|) 
to be rewritten in recursive form 

= P(+)* i (+)+F(-)* i (-), (15) 
= F(+)*4+)+F(-)* i (-), (16) 

either of which allows one to construct the full solution 
for the ideal wire, once the boundary values are set, cf. 
Eq. ©. 



C. Scattering problem 

In a non-ideal quantum wire the potential is not peri- 
odic, which means that we have to solve the Schrodinger 
equation of Eq. 10 with depending upon the posi- 
tion i along the wire. The non-ideal region (the scatter- 
ing region) is supposed to be finite, spanning the cells 
i = 1, . . . , Sm^ The left and right leads are ideal wires, 
spanning the cells i = — oo, . . . , and i = S + 1, . . . , oo, 
respectively. In the ideal wires does not depend on 
the position of the cell. However, the left lead can be 
different from the right one, so we use the subscript L(R) 
to denote the former (latter), i.e. Hi = Hl, i < 1 and 
Hi = Hr, i > S. A schematic picture of the structure 
is shown in Fig. EJa). We solve Eq. 10 over the whole 
space, i = — oo, oo, making use of the ideal wire solu- 
tions of the previous section to reduce the problem to 
essentially the scattering region only, see Fig.^Jb). 

For the solution in the left lead the recursion relation 
of Eq. (|To|> can be used. This gives for the cell i = — 1 

¥_i = F l (+)* (+)+Fl(-)*o(-) 

= [F l (+)-F l (-)]*o(+)+Pl(-)*o, (17) 



using Eq. 0. The vector ^o(+) describes a wave coming 
in from the left. In a scattering problem this vector fixes 
the boundary condition. Eq. I|17fl allows Eq. JSJ for i = 
to be written as 



where 



(SI-H )*o + Bt* 1 =Q* (+), (18) 

Ho = H L -BF L (-) 

Q = B[F l (-)-F l (+)]. (19) 



For the solution in the right lead we use the recursion 
relation of Eq. (|15fl , which gives for the cell i = S + 2 



* S+2 = F R (+)* S+ i(+). 



(20) 



Here we have assumed that in the right lead we have only 
a right going wave, which corresponds to the transmitted 
wave. Eq. (|20|l allows Eq. |J5J| for i = S + 1 to be written 

as 



where 



(£I-H s+ i)* s+ i +B*. 



Hs+i — Hr — B^Fr(+). 



0. 



(21) 



(22) 



Eqs. I|19f) an d (|22ll take care of the coupling of the 
scattering region to the left and right leads. Eq. (Jsj for 
i = l,...,S plus Eqs. (fTHfl and form a complete set 
of equations from which the vectors i = 0, . . . , S + 1 
can be determined describing the waves in the scattering 
region. 

The scattering reflection and transmission coefficients 
can be deduced from the amplitudes immediately left and 
right of the scattering region, i.e. \&o and &s+i- If we let 
the incoming wave consist of one specific mode, 1 4 , o(+) 



ul. 



i.e. a, 



.(+) 



in Eq. 10, then the generalized 



reflection and transmission probability amplitudes r„<„ 
and t n >„ are defined by 

M L 

*o(-) = ^2 UL,n'(-)»Vn 
n'=l 

Mr 

*s+i(+) = J2 ur <"' (+)*«'»• ( 23 ) 

n' = l 

Note that at this stage we include all evanescent and 
propagating modes since these form a complete set to 
represent the states in the leads. We assume the lead 
states to be amplitude normalized. 

The reflection and transmission probability amplitudes 
r n i n and t n i n between all possible modes form a Ml x 
Ml matrix R and a Mr x Ml matrix T, respectively. 
All elements of these matrices can be found in one go 
by defining a N ts x Ml matrix of all possible incoming 
modes, i.e. 



C (+) = U L (- 



(24) 
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Analogous to Eq. 112,'Sll one then has 

Co(-) = Co-C (+)=U L (-)R 
C s+ i(+) = C 5+1 =U R (+)T. (25) 

Eqs. (jTHjl . © and l(2*TJl then become 

(M-H )C + BtC 1 = QU L (+) 
(M-HiJCi+BCi-i+BtCi+i =0 (26) 
( J BI-H s+1 )C s+1 +BC s = 0, 

i = 1, . . . , S. Solving this set of equations for d; i = 
0, . . . , S + 1 gives all possible waves. From Eq. lj2*5l) one 
can then extract the generalized reflection and transmis- 
sion matrices R and T. An efficient technique for solving 
the equations is discussed in appendix IA 41 

In order to calculate the total transmission one has 
to select the transmission matrix elements that refer to 
propagating modes and discard the ones that refer to 
evanescent modes. This is easy, since the propagating 
modes have |A| = 1, see the discussion above Eq. JSJ. 
The total transmission of Eq. is then given by 

m L ,mR 

T (E)= V ^\t n 'n\\ (27) 
n— l,n — 1 

where VR,. n > and are the velocities in the x-direction 
of the right propagating waves in the right and left lead 
in the modes n' and n, respectively, and tol, tir are the 
number of such modes. Introducing the velocities results 
from flux normalizing the modes, which is required by 
current conservation^ The velocities are given by the 
expression 

v n = -^Im [A„ut Bt u „] , (28) 

where subscripts L(R) need to be added for the left(right) 
leads. Eq. (|28|1 is derived in appendix IA 51 The sign of 
the calculated velocities is used to distinguish right from 
left propagating modes. 

D. Computational costs 

Our computational method can be summarized as 
follows. First, Eq. Q is solved in its reduced form, 
Eq. HA8jl . to obtain the modes for both leads. The 
computing costs of this step scale as Nh, where 
N id = max(2A,i — N) ■ W y ■ W z , see appendix lA"2l 
These costs are small compared to the costs of solving the 
scattering problem. The next step involves the selection 
of the physically relevant modes u m and separate them 
into left (+) and right (— ) going modes. The velocities 
are calculated using Eq. (|2*5)l and are used to distinguish 
left from right propagating states. Evanescent states are 
classified as growing (+) or decaying (— ) on account of 
their eigenvalue. Subsequently, the F-matrices are con- 
structed, Eq. (|13|) . and the matrix elements that define 
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FIG. 3: (Color online) The potential of Eq. |J2§) in the air- 
plane for the cases (a) Vb = and (b) Vo = Vi . 

the boundary conditions on the scattering region are set, 
see Eqs. (|19f) and l|22|) . The computing costs of these 
steps are minor. 

The transmission matrix T is obtained by solving 
Eq. (|26|l using the algorithm of appendix IA 41 This is 
the most time consuming step. It scales as S ■ N? s , where 
N IS = L ■ W y ■ W z is the number of grid points in the unit 
cell and S is the number of unit cells in the scattering 
region. Note that the scaling is linear with respect to 
the size S of the scattering region, which means that this 
algorithm can be classified as O(N). Finally, the total 
transmission and the conductance can be obtained from 
Eqs. (|77jl and ffl. 



III. RESULTS 

A. Numerical tests 

In order to test the accuracy of our method we consider 
a system described by the model potential 

V(r) = Vq [cos(27re/a) + cos(27ry/a) + cos(27rz/a)] 



cosn {wx/a) 

The Vq term describes an ideal wire by a simple three 
dimensional periodic potential with periods a x = a y = 
a z = a. The V\ term describes a barrier in the propaga- 
tion direction and is a simple model for a tunnel junction. 
The potential is plotted in Fig. [31 We solve the scattering 
problem for this system numerically in three dimensions 
by the method outlined in Sec. [H] Our results can be 
verified, however, since this potential is in fact separable 
and limiting cases can be solved analytically. The solu- 
tions in the y- and z-directions are Mathieu functions^ 
If V\ = then the solutions in the x-direction are also 
Mathieu functions. If V\ ^ but Vb = the scattering 
problem can be solved analytically^ Finally, if V\ ^ 
and Vb ^ the solution in the x direction can be obtained 
using the separability of the potential and a standard nu- 
merical solver for the resulting ordinary differential equa- 
tion in the a;-directioni2& In the following the latter will 
be called the "exact" numerical solution. 
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FIG. 4: The calculated wave number k x (in units of n/a) 
as function of the energy E (in units of Vb) for an ideal 
wire. The points indicate the numerical results obtained with 
L, W y , W x = 8 and N = 4. The solid line indicates the exact 
solution of the Mathieu problem. 

As a first test we consider an ideal wire, i.e. V\ = 
in Eq. I|29|) . The potential is separable and we can write 
the energy as E{k x ,k y ,k z ) = e nx {k x ) + e ny (k y ) + e nz (k z ), 
where e„(fc); k = —tt/cl, . . . ,7r/a; n = 0, 1, . . . are the 
eigenvalues of the Mathieu problem^ Fig. 0] shows 
part of the analytical band structure for {k x , k y , k z ) — 
(0, 0, 0) — ► {it /a, 0, 0). It essentially consists of a superpo- 
sition of one-dimensional band structures t nx {k x ) offset 
by energies e ny {0) + e n2! (0); n y ,n z = 0, 1, . . .. 

The numerical band structure is obtained by solving 
Eq. Q in its reduced form, Eq. (|A8|) . To obtain the re- 
sults shown in Fig. 0] we set k|| = (k y ,k z ) = (0,0), cf. 
Eq. (IA4I) . and determine the eigenvalues A in Eq. J7J) as 
a function of E. For the propagating states one can write 
A = exp(ik x a). Plotting the calculated wave number k x 
as function of the energy E then allows to compare the 
results with the analytical band structure. The numer- 
ical results shown in Fig. 0] are obtained using a grid of 
L, W y , W z = 8 points per period and a FDA with N = 4. 
Although this grid is relatively coarse, we obtain a rela- 
tive accuracy on k x of 10~ 3 . 

This perhaps surprising accuracy is entirely due to the 
use of a high order FDA. To illustrate this, Table[I]shows 
the convergence of k x at a number of energies E as a func- 
tion of the order N of the FDA and the number of grid 
points L. These particular results were obtained using 
the separability of the potential and solving the problem 
numerically in the x-direction only, while using analyti- 
cal solutions for the y- and z-directions. For N = 6 and 
L = 14 the results are converged to within 10 -7 of the 
exact result. This is in sharp contrast to the results ob- 
tained with a simple first order (N — 1) FDA, where a 
similar convergence can only be obtained at the cost of 
using two orders of magnitude more grid points. Using 
such a large number of grid points in three dimensions 
is entirely prohibitive because of the high computational 



TABLE I: k x (E) (in units of 7r/a) at values of E (in units 
of Vb) in the lowest two bands and in the first band gap of 
Fig-El m the band gap we find k x = 1 + lK x - 



N 


L 


M-0.6) 


Ml.0) ^(0.3) 


1 


7 


0.694906 


0.283091 0.283549 




10 


0.608958 


0.322712 0.304238 




14 


0.571387 


0.342011 0.312259 




100 


0.533341 


0.359187 0.319658 




1000 


0.533084 


0.359425 0.319688 


4 


5 


0.544347 


0.355049 0.317777 




7 


0.533962 


0.358927 0.319476 




10 


0.533149 


0.359389 0.319673 




14 


0.533087 


0.359425 0.319687 


6 


7 


0.533228 


0.359337 0.319658 




10 


0.533086 


0.359425 0.319688 




14 


0.533082 


0.359428 0.319688 


Exact 




0.533082 


0.359428 0.319688 



costs involved. For example, aiming at a moderate accu- 
racy of 10~ 2 , it is observed that for TV = 4 and L = 5 the 
results are markedly better than for N = 1 and L = 14. 
Yet in a three dimensional calculation, without using the 
separability of the potential, the computing time required 
for the latter is two orders of magnitude larger than for 
the former. It means that in order to solve a general non- 
separable three dimensional problem with reasonable ac- 
curacy and computational costs, it is vital to use a high 
order FDA. 

Next we consider the scattering problem and calculate 
the total transmission for the case where V\ ^ 0. The 
size of the scattering region is set to Sa and outside this 
region the scattering potential (the last term of Eq. (|29(l ^l 
is set to zero. With S = 6 the results are extremely well 
converged. As an example we have calculated the trans- 
mission at normal incidence, i.e. k|| = (k y ,k z ) — (0,0). 
The crosses marked Vq = in Fig. [S] represent the nu- 
merical results for the transmission of the correspond- 
ing potential, obtained with a L, W y , W z — 8 grid and a 
N = 4 FDA. This scattering problem can also be solved 
analytically, 37 and the analytical and numerical trans- 
mission probabilities agree within 10~ 4 . 

Fig. [3] also shows the transmission for the case where 
Vq = V\ as calculated numerically using the same param- 
eters as before, i.e. S = 6;L,W y ,W z = 8; N = 4. This 
scattering problem can only be solved semi-analytically; 
Mathieu solutions are used in y- and z-directions, and 
the (ordinary) differential equation for the x-direction is 
solved "exactly" using an accurate standard numerical 
solver* 3 -^ Again the "exact" and numerical transmission 
probabilities agree within 10 -4 . Compared to the Vb = 
case it is observed that the influence of the periodic po- 
tential of the leads upon the transmission is large. For 
Vq = V\ the electronic states in the leads are far from free 
electron-like, see Fig. In particular, the transmission 
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FIG. 5: (Color online) The total transmission as function of 
energy (in units of Vi) for the two cases Vb = (crosses) and 
Vb = Vi (dots). In both cases kii = (0,0). The solid lines 
represent the analytical solution for Vb = and the "exact" 
numerical solution for Vb =Vi. 




FIG. 6: (Color online) The total transmission as function of 
energy (in units of Vi) for the two cases kn = (0,0) (dots) 
and k|| = (0.47, 0.21)n/a (triangles). In both cases Vb = Vi. 
The solid lines represent the "exact" numerical solution. 
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0.1321 
0.1319 
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FIG. 7: The total transmission T as function of the grid size 
L = W y — W z for a simple N = 1 FDA,'+' top curve, and 
for a N = 4 FDA, 'x' bottom curve. The horizontal line 
represents the "exact" value T = 0.132. The inset shows the 
N — 4 curve on a finer scale for T. 



and a high order N = 4 FDA. The results shown are for 
one particular k|| = (k y ,k z ) — (0.47, 0.21)7r/a and energy 
E = 0.895 Vb, but the convergence at other k||-points 
and energies is very similar. The number of propagating 
channels at this k|| -point and energy is two, but the total 
transmission is only T — 0.132, which means that the 
barrier is largely reflecting. We conclude that the accu- 
racy of the three-dimensional calculation depends very 
strongly upon the order of the FDA. For N = 1, L = 15, 
the transmission is converged on a scale of 10 -2 only, but 
for N = 4 it is converged on a scale of 10~ 3 already for 
L — 8, see the inset of Fig. El A high order FDA thus en- 
ables the use of a much coarser real space grid. Since the 
computational costs scale with the number of real space 
grid points N IS — L 3 as N^ s = L 9 , this demonstrates the 
strength of using a high order FDA. 



B. Sodium atomic wires 



drops to zero if the energy is inside a band gap, because 
there are no lead states of that energy. 

The numerical calculations accurately capture the 
transmission curve over a large energy range, as is shown 
in Fig. |SJ The transmission generally increases with 
energy due to the increasing number of channels, sec 
Fig. 0] Since the density of states peaks at the band 
edges, the transmission peaks at the corresponding en- 
ergies. The transmission depends very much upon kp as 
can be observed in Fig. where the transmission for nor- 
mal incidence, k|i = (0,0), can be compared to that for 
k|| = (0.47, 0.21)7r/a (an arbitrary point in the Brillouin 
zone). The difference between the two curves can be eas- 
ily understood from the band structure of the leads. In 
particular, for k|| = (0.47, 0.21)7r/a there are no band 
gaps for E > 0.14Vb. 

To demonstrate the convergence of the numerical cal- 
culations, Fig.|2shows the total transmission as function 
of the sampling density L = W y = W z for a simple N = 1 



We have calculated the electronic transport in sodium 
atomic wires as examples of more complex systems. Our 
model of a sodium wire consists of left and right leads 
composed of bulk (bcc) sodium metal terminated by a 
(100) surface, connected by a straight wire of sodium 
atoms, as is shown in Fig- EI The atoms in the leads are 
positioned according to the bcc structure of bulk sodium, 
with the cell parameter fixed at the experimental value of 
7.984 aom& The atoms in the wire are fixed at their (bulk) 
nearest neighbor distance of 6.915 ciq. Since geometry 
relaxation at the Na(100) surface is very smallj^H and 
calculations using jellium electrodes have shown that the 
conductance of a sodium wire is not very sensitive to 
its geometry^i we have refrained from optimizing the 
geometry. Perpendicular to the wire we apply periodic 
boundary conditions using a 2 x 2 lateral supercell, which 
has a lattice parameter of 15.968 ao- 

If DFT is used to model the electronic structure, Eq. 
@ corresponds to the Kohn-Sham equation. The one- 




FIG. 8: (Color online) Structure of an atomic wire consisting 
of 4 sodium atoms between two sodium leads terminated by 
(100) surfaces. The scattering region is bounded by the verti- 
cal lines and the lateral supercell by the horizontal lines. Bulk 
atoms are indicated by yellow (light grey) balls and atoms in 
the scattering region by blue (dark grey) balls, respectively. 



electron potential V(r) in this equation is then given 
by the sum of the nuclear Coulomb potentials or pseu- 
dopotentials, and the electronic Hartree and exchange- 
correlation potentials. The latter two depend upon the 
electronic charge density. In linear response the charge 
density remains that of the ground state, allowing the 
electronic potentials to be obtained from a self-consistent 
ground state calculation. In these calculations we employ 
supercells containing a slab of 13 layers to represent the 
bulk and surface of the leads, and a wire of n atoms, see 
Fig. [5] We use the local density approximation^ and 
represent the ion cores of the sodium atoms by a local 
pseudopotential4*i The valence electronic wave functions 
are expanded in a plane wave basis set with a kinetic en- 
ergy cutoff of 16 Ry. The lateral Brillouin zone is sampled 
with a 8 x 8 k|| -point grid, using a temperature broaden- 
ing with kT el =0.1 eVii 

We want to calculate the conductance for various 
lengths of the atomic wire, so for each length n we per- 
form a self-consistent supercell calculation to generate 
the one-electron potential. By expressing the latter in 
a plane wave basis, Fourier interpolation can be used to 
obtain a representation on any real space grid required 
for the transport calculations, cf. Eq. (J3J. An example 
of the effective potential for valence electrons of a sodium 
atomic wire is shown in Fig. [5] 

For the transport calculations we use a scattering re- 
gion comprising the atomic wire and the surface regions 
of the left and right leads. Both surface regions consist of 
5 atomic layers, see Fig.|SJ Periodic boundary conditions 
and a 2 x 2 lateral supercell perpendicular to the wire 
are applied. The potential in the scattering region is ex- 
tracted from the slab calculations. Outside the scattering 
region we assume that the leads consist of bulk sodium. 
The potential for the leads and the value of the Fermi 
energy are extracted from a bulk calculation. The aver- 
age bulk potential is lined up with the average potential 
in the middle of the supercell slab4^ We have checked 
that the spatial dependence of the bulk potential is vir- 
tually identical to that of the potential in the middle of 
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FIG. 9: (Color online) The effective potential for valence elec- 
trons in the xy-plane of the sodium atomic wire shown in 
Fig |S| Most prominent are the strongly repulsive core re- 
gions and the attractive valence regions of the atoms. The 
difference between the maximum and minimum values of this 
potential is 32.4 eV. The Fermi level is at 8.5 eV above the 
potential minimum. 



the slab. This means that the connection between the 
leads and the scattering region is smooth. There are no 
discontinuities in the potential in the boundary regions 
that could cause spurious reflections. 

Fig.llOlshows the calculated conductance G(Ep) at the 
Fermi level as a function of the length n of the sodium 
atomic wire. The conductance is calculated using a grid 
spacing along x,y and z directions of h Xty<z — 0.67 ao 
(giving L = 6, W y = W z = 24) and a FDA of order 
N = 44& The 2x2 lateral Brillouin zone is sampled by 
a 12 x 12 uniform k||-point grid. With these parame- 
ters the calculated conductances have converged to well 
within 10 -3 Go (where Go = e 2 / (7rS))4£ All wires have a 
conductance close to unity, except the one-atom (n = 1) 
wire. The high conductance for the one-atom wire is fore- 
most due to tunneling between the left and right leads 
through vacuum. The latter can be calculated by omit- 
ting the wire and otherwise keeping the geometry fixed. 
Tunneling through vacuum leads to a conductance of 2.20 
Go per 2x2 surface cell. Vacuum tunneling decreases fast 
with the distance between left and right lead surfaces. At 
a distance corresponding to the two-atom wire it gives a 
conductance of 0.047 Go; at distances corresponding to 
longer wires this conductance is negligible. Subtracting 
these vacuum tunneling values gives the lower curve in 
Fig.Uni 

With the exception of the one-atomic wire, the conduc- 
tances are close to unity. This is perhaps not surprising, 
since within a tight-binding model atomic wires consist- 
ing of a monovalent atom like sodium are expected to 
have one open channel. For perfectly transmitting con- 
tacts this would give a conductance of 1 Go- 48 Our cal- 
culated conductances for n > 1 are less than 15% smaller 
than this value, demonstrating that this transmission is 
indeed very high. The conductance of the one-atomic 
wire, relative to the vacuum tunneling conductance at 
this distance, is significantly lower than unity. The elec- 
tronic structure of a single atom between two electrodes is 
substantially distorted from the simple single open chan- 
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FIG. 10: Top curve: conductance G{Ef) at the Fermi level 
(in units of e 2 /irh) of a sodium atomic wire as a function of 
the number of atoms n in the wire. Bottom curve: conduc- 
tance of a sodium atomic wire relative to vacuum tunneling 
conductance between two electrodes without wire. 



nel models 

On a finer scale we find evidence of an even-odd 
oscillation in the conductance obtained in previous 
studiesi£i22i£2i£iiS The conductance for wires with n even 
tends to be lower than for those with n odd. In simple 
tight-binding terms odd-numbered atomic wires have a 
non-bonding level that tends to line up with the Fermi 
level of the leads, which gives a high transmission. For 
an even-numbered atomic wire on the other hand the 
Fermi level tends to fall in the gap between the bonding 
and anti-bonding levels of the wire, resulting in a lower 
transmission^ The size of the even-odd oscillation in the 
conductance depends of course upon the nature of the 
contacts between the atomic wire and the leads. Good 
contacts broaden the levels of the atomic wire into wide 
resonances, which tends to suppress the even-odd oscil- 
lation. Our wires have good contacts, but the even-odd 
oscillation remains distinctly visible. 54 



IV. SUMMARY 

We have formulated and implemented a new numerical 
technique for calculating electronic transport in quan- 
tum wires and tunnel junctions in the linear response 
regime, starting from Landauer's scattering formalism. 
It is based upon a real space grid representation of the 
scattering problem. Dividing space into left and right 
ideal leads and a scattering region, the problem is solved 
by wave function matching (WFM). First all propagat- 
ing and evanescent Bloch modes of the leads are calcu- 
lated. Subsequently the states in the scattering region 
are forced to match to the Bloch modes of the leads. 
This directly leads to the transmission matrix, which con- 
tains the transmission probability amplitudes between all 
modes of the left and right leads, and to the conductance. 
The computing costs of this algorithm scale linearly with 
the size of the scattering region. 



It is shown that the use of a high order finite dif- 
ference approximation for the kinetic energy operator 
leads to a high accuracy and efficiency This is demon- 
strated for a model potential by benchmarking the tech- 
nique against analytical and numerically "exact" solu- 
tions. The method is then applied to calculate the con- 
ductance in sodium atomic wires, where the potential in 
the wire and in the bulk sodium leads is obtained from 
self-consistent DFT calculations. 
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APPENDIX A: COMPUTATIONAL DETAILS 

1. H and B matrices 

In this section the matrices and B, introduced in 
Sec. Ill Al are presented explicitly. In order not to compli- 
cate the notation the subscript i is dropped; all quantities 
refer to a single cell. For a wire with a finite cross-section 
in the yz plane, the matrix H is real and symmetric and 
has the form 



/ hi -fr 

-Pi h 2 



H 



-Pn 
~Pn-i 





-Pn 



\ 




/ 



... ... h L _i -Pi 
V ... ... —Pi h L 

(Al) 

Here N is the order of the finite difference formula used, 
see Eq. ©. We assume that the x-axis is in the direction 
of the wire. L is the number of grid points in the x- 
direction of the unit cell defined by the periodic potential. 

The submatrices h„ and p n are of dimension W y x W z , 
which is the number of grid points in the cross sec- 
tion of the wire. Denoting (k,l) = k + (I — l)W y ;k 

I W v ; I — 1, . . . , W z as the compound index covering 

the grid points in the cross section, the non-zero elements 
of these matrices are easily derived from Eq. J1J . 



(h,-)(M),(fc,0 
(hj)(M),(fc+n,0 
(hj)(&,;),(W+n') 

(l3j)(k,i),(k,i) 



V*,W-(*o+*8 + *o) 
-tl , n ± 



(A2) 
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where -N < n,n' < N and 1 < k, k + n < W y ;l < 
1,1 + n' < W z ;l < j < N. Note that in writing down 
these matrices we have assumed that N < L, W y , W z . In 
practical calculations on realistic systems this will always 
be the case. 

The matrix B has the same dimension as H, but it is 
upper triangular 







B 



(3 N (3n-i 
p N 





Pi \ 

P2 



P 



N 



\0 ... 



(A3) 



o / 



For a wire that is periodic in the yz plane, the wave 
functions in Eq. Hi must obey Bloch conditions. That 
is, ®j,k+W v ,l = e^ a "% fe>! and tyj, k ,l+w. = e ifc * a **i, fc ,;, 
where a y , a z are the periods in the y and z directions, and 
{k y ,k z ) = k|| is the Bloch wave vector in the yz plane. 
These Bloch conditions in the yz plane can be taken into 
account by defining the blocks 

{^'j){kj),{k+w v +n,i) = -tleT ih * a \ n = -N,...,-k 

{h'j)(k,l),(k-W y +n,l) = -t\^ a \ n = W y -k,...,N 

0^'j)(k,l),(k,l+w z +n') = -t^,e~' lkzaz , n' = —N, . . . , —I 

(h' J 0(fc,O,(*,i-w,+n') = -t z n >e ikza % n' = W z -l,...,N. 

(A4) 

The matrix H(ky), which is obtained by substituting hj 
by hj ;+ h'j ; ; j — 1, . . . , L in Eq. (|A1|) describes a wire that 
is periodic in the yz plane with solutions corresponding to 
a Bloch vector k|| . This matrix is (complex) Hermitian. 



2. Ideal wire 

For an ideal wire, which has a periodic potential along 
the wire, Eq. Q has to be solved to find the propagat- 
ing and the evanescent waves. The precise form of the 
submatrices in Eqs. (|A1|) and (|A3|) is not important in 
the following discussion. For ease of notation we only 
mention the dimensions L (the number of grid points in 
the x-direction) and N (the order of the finite difference 
expression) explicitly and treat the wire as quasi one- 
dimensional. To find the dimensions of the matrices in 
the three-dimensional case, one simply has to multiply 
the dimensions mentioned below by W y x W z . 

Eq. Q is a generalized eigenvalue problem of dimen- 
sion 2L. Because the matrix B is singular it has a number 
of trivial solutions A = and A = oo. By using a parti- 
tioning technique we will eliminate these trivial solutions 
and reduce the problem to the 2N non-trivial solutions. 
The key point is to split the vectors into two parts 
containing the first L — N and last N elements, respec- 
tively. The two parts are denoted by the subscripts 1 and 



2. Splitting the matrices H and B in the same way one 
gets 




Hn H12 
H21 H22 



; B = 



B 12 \ 
B 22 J ' 

(A5) 

Note the special form of the matrix B. 

This splitting allows Eq. JJJ to be written in the form 



AB 21 



/ EIu — Hn 

— H21 

In 



*M \ 

V *i-l,2 J 



Eh 



— H12 
H22 -f 


I22 



AB 22 






B12 \ 





B22 


-AIn 








-AI22 / 



= 0, 



(A6) 



From this expression it is clear that the component 
only enters the problem in a trivial way as 
^i-1,1 = 1/A x Hfjji. It can be eliminated by deleting 
the third row and column in Eq. (|A6p . 

Furthermore the first row of the matrix does not de- 
pend upon the eigenvalue A. Writing out the multiplica- 
tion for the first row explicitly, one finds an expression 
for 

* M = (EI n - Hn)" 1 (H^*.^ - B^-1,2) . (A7) 



This can be used to eliminate ^^1 from Eq. I|A6|I to 
arrive at the equation 



- A 





with 



An — EI22 — H22 — H21 (Sin — Hn) 1 H12 
A12 = H 2 i (Sin — Hn) 1 B12 



Sn 
S12 



'B 2 2 — B 21 (Sin — Hn) 1 H12 
"B 2 i (Sin — Hn) 1 Bi 2 . 



(A9) 



Eq. I|A8(1 is a generalized eigenvalue problem of di- 
mension 2N that can be solved using standard numer- 
ical techniques^ In general it gives 2N eigenvalues A m 
and eigenvectors u m . As mentioned in the text, some 
of these solutions are non-physical ^iSM others represent 
extremely fast growing or decaying waves. Both of these 
classes of unwanted solutions are easily filtered out by 
demanding that 1/8 < |A| < S, where S is some thresh- 
old value. We use this criterion to select the physically 
relevant solutions, which are then separated into M right- 
going and M left-going solutions. These are used to con- 
struct the matrices of Eqs. © and l|ll|) which contain all 
the information required to describe the ideal wire. 
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The computational cost of solving Eq. (|A8|) scales as 
(27V) 3 , whereas the cost of computing the matrices of 
Eq. (|A9(I basically scales as (L — N) 3 (which is the cost of 
the matrix inversion involved) . Depending on the relative 
sizes of L and N one of these two steps is dominant. 



making use of Eqs. (|11|) and i|A10|) - (|A15|) . In a similar 
way a solution to Eq. (|14|l is formed by 

F = UA _1 U 

M 

= |u m )A m 1 (a m |. (A18) 



3. F matrices 

In this section explicit expression for the matrices F 
and F are given, see Eqs. <|13[) and l|14|) . Following 
Eq. JSJ), we denote the propagating and evanescent modes 
of the ideal wire by u m ; m — 1, . . . , M, where M < N Ta 
and N IS is the dimension of the vectors. For clarity of 
notation we omit the labels ± for right- and left-going 
modes here. As in Eq. JHJ we form the N rs x M matrix 



U = (m ••■ u M ) 

Uu ... Uim 



(A10) 
(All) 



UN r! ,M 



The mode vectors u m are in general non-orthogonal and 
we can form the M x M (positive definite) overlap matrix 
with elements 



Srnn — U m lln — (u m jll n ). 

This allows us to construct the dual basis u r , 
1,...,M 



(A12) 



M 



»=1 



with properties 



(u m |u„) = (u m |u„) 
Now define the M x iV rs matrix 
U = (ur-UMjt 



(A13) 



(A14) 



u 



li 



'1M 



N tB l 



l N IB M 



(A15) 



U is called the pseudo- inverse of U; note that UU = Im, 
where Im is the M x M identity matrix^ 
Defining the matrix 



UAU, 



(A16) 



it is easy to show that it is a solution to Eq. (|13fl . F is in 
fact a matrix that projects onto the space spanned by the 
modes, as is easily demonstrated by writing Eq. I|A16I) as 



M 

E 

m— 1 



|u m )A 



(A17) 



Note that F = F" 1 
N x W x x W y < N n 
never be the case. 



only if M = N rs , but since M < 
(see the previous section) this will 



4. Scattering problem 

The scattering problem is described by Eq. Ij26(l . It is 
conveniently written in matrix form as 



/A B' 
B Ai Bt 
B A 2 



V o 

with 



Co 
C 2 








. . . A s+ i / V c s+i / 



(A19) 



A 
A, 

As+i 
C 
Cs+i 
D 



Ei-n 

EI U t , i = l,... 
EI - H s+ i 
U L (+) + U L (-)R 
U r (+)T 
QU L (+). 



s 



(A20) 



All the blocks A-D are A^ rs x N rB matrices. Eq. (|A19|) 
represents a set of linear equations, which can be solved 
directly using a standard algorithm. However, the di- 
mension of this problem is N to t = N Ta ■ (5 + 2), which 
can be rather large. Since the computing cost scales as 
N 3 ot the direct route is not very practical. 

It is however quite straightforward to construct an al- 
gorithm for which the computing cost scales as N 3 S -5, i.e. 
only linearly with the size 5 of the scattering region. One 
has to make optimal use of the block tridiagonal form of 
the matrix in Eq. (|A19|) . The algorithm is a block form 
of Gaussian elimination. The first (and most time con- 
suming) step of this algorithm is transforming the matrix 
into upper block triangular form by iteration 



A' = A ; D' =D; 

A'i = Ai- BA'T^Bt 
D'i - -BA'^D'j-i 



i = 1,...,5+1.(A21) 



The inverse matrices A'j_j in this algorithm are actually 
not needed explicitly. Instead at each step one solves the 
sets of linear equations 



A'i_iBi = Bt ; A',_iD, = D', 



(A22) 
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by a standard algorithm, i.e. LU decomposition of A'j_i 
followed by back substitution, to obtain the matrices B, 
and D,;. 55 This allows the steps in Eq. (|A21|) to be rewrit- 
ten as 

A'j - A, - BBj ; D' t = BD, (A23) 

The solution to Eq. l)A19j) can now be found by back 
substitution 

Cs+i = A' g+1 D's +1 

Ci = D i+ i - B 4+1 C l+1 } i = S, . . . , 0. (A24) 

Again one does not need A'T^ explicitly, but like 
Eq. (|A22|> one can solve the equivalent set of linear equa- 
tions. The reflection and transmission matrices R and T 
can be extracted using the special form of the matrices 
C and C s +i, see Eq. (|A20j) . 

Very often one is only interested in the transmission 
matrix. In that case one only uses the first step of the 
back substitution, Eq. (|A24I) . which can be written as 

A' S+1 U R (+)T = D' S+1 . (A25) 

This is a set of linear equations for the transmission prob- 
ability amplitudes T, which can be solved using a stan- 
dard numerical techniques S 

The time consuming steps consist of solving Eq. I|A22J| . 
the computing costs of which scale as N^Ml Using 
Eq. (|A"23)) in Eq. fMQ requires performing S + 1 of such 
steps and subsequently solving Eq. ljA25(l scales as iV r 3 s . 
Note that the full algorithm scales linearly with the size 
S of the scattering region. 



This quadratic eigenvalue equation of dimension N rs is 
completely equivalent to the linear problem of dimension 
2N rs of Eq. 0. If u m is a right eigenvector of Eq. (JA26|) 
belonging to the eigenvalue A m , then by complex conju- 
gation of this equation one shows that is a left eigen- 
vector belonging to the eigenvalue 1/A^. For a prop- 
agating state |A m | = 1, so A m = 1 / , which means 
that these left and right eigenvectors belong to the same 
eigenvalue. 

We now start from 



A m u]; i ( J EI-H)u m + u] n Bu m + A^uJ ri B t u m = 0, (A27) 



and take the derivative d/dE of this expression. All the 
terms with du m /dE and du^/dE drop out, because u m 
and obey Eq. I|A26|) and its complex conjugate, re- 
spectively. The remaining terms can be collected and 
slightly rewritten using Eq. i|A26[l : the result is 

(A^uJ^Bu™ - A^u^B+u™) + A m u,|„u r „ = 
_ 2l ^il m (A TO u+ 1 B t u ro ) + A m = 0, (A28) 

where the last line is obtained by making use of A" 1 = A* n 
and the fact that the vectors are normalized, uj„u ra = 1. 
Eq. (|A28() yields an expression for dX m /dE. For propa- 
gating states A m = e lfcx<Ix and thus 



5. Velocities 

In this section we give a short derivation of the expres- 
sion for the velocities, Eq. i|28|) . It is straightforward to 
show that the vectors u,„ of Eq. © are a solution of the 
quadratic eigenvalue equation 

X m {EI - H)u m + Bu m + A, 2 ri Bt Um - 0. (A26) 



dE ia x X m dE 



The usual definition of the Bloch velocity v n = 
H~ 1 dE / dk x and the expression for d\ m /dE extracted 
from Eq. (|A28|I then give the expression for the veloc- 
ity of Eq. 
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